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Rapidly computable viscous friction and no-slip 

rigid contact models 

Evan Drumwright 


Abstract —This article presents computationally efficient al¬ 
gorithms for modeling two special cases of rigid contact— 
contact with only viscous friction and contact without slip—that 
have particularly useful applications in robotic locomotion and 
grasping. Modeling rigid contact with Coulomh friction generally 
exhibits 0{n^) expected time complexity in the number of contact 
points and 2°^"^ worst-case complexity. The special cases we 
consider exhibit 0{m^ + m^n) time complexity (m is the number 
of independent coordinates in the multi rigid body system) in 
the expected case and polynomial complexity in the worst case; 
thus, asymptotic complexity is no longer driven by number of 
contact points (which is conceivably limitless) but instead is 
more dependent on the number of bodies in the system (which 
is often fixed). These special cases also require considerably 
fewer constrained nonlinear optimization variables thus yielding 
substantial improvements in running time. Finally, these special 
cases also afford one other advantage: the nonlinear optimization 
problems are numerically easier to solve. 

I. Introduction 

Dynamic robotic simulation; grasp planning; and, increas¬ 
ingly, locomotion planning and control employ rigid contact 
models with dry (typically Coulomb) and wet (viscous) fric¬ 
tion. These contact models yield an effective tradeoff between 
computation speed and physical accuracy. While rigid contact 
models are far faster than, e.g., elastodynamic finite element 
analysis, they still require heavy computation: the expected 
time complexity for such models is O(n^) in the number of 
contact points. Additionally, the number of contact points input 
to the model is conceivably limitless. This issue is not just 
theoretically interesting: Wang Il25l reports that solving the 
contact problem absorbs up to 90% of computation time when 
simulating a scenario for the DARPA Robotics Challenge 
using ODE m. 

Roboticists are often content to use rigid contact models 
without Coulomb friction for computational expediency. For 
example, one may wish to model locomotion or effect sim¬ 
ulated grasping without observing slip; roboticists studying 
legged locomotion often predicate their models on no slip 
occurring, for example. If slip is desirable, purely viscous 
friction might yield a suitable model if, for example, a robot is 
walking on a wet surface. This article presents computationally 
efficient methods for both of these special cases. 

These special cases provide the following computational and 
modeling advantages: (i) time complexity goes from worst- 
case exponential (the worst-case complexity of solving rigid 
contact problems with Coulomb friction Gull] using Lemke’s 
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Algorithm ffSl ) to worst-case polynomial in the number of 
contact points; ( 2 ) significant reduction in the number of 
nonlinear optimization problem variables; and ( 3 ) a positive- 
semi-definite-matrix linear complementarity problem (LCP), 
in place of a copositive-plus LCP, which is demonstrably 
easier to solve m (i.e., the solver is less likely to fail due 
to numerical errors) and permits the use of general algorithms 
for solving convex optimization problems. 

Finally, we provide an algorithm that yields 0{m^ + m?n) 
expected asymptotic time complexity on these two contact 
models, where m is the number of independent coordinates in 
the multi rigid body system. This algorithm therefore provides 
a means to make complexity more dependent on the number 
of independent coordinates in the system (this number remains 
constant except in the unusual case in which bodies are 
inserted into the simulation) than on the number of contact 
points (which is conceivably unlimited). 

II. LCPS, NCPS, AND MLCPS 

A LCP, or linear complementarity problem, (r, Q) signifies 
the problem: 

w = Qz -f r 
w > 0 
z > 0 
z^w = 0 

for unknown vectors z,w € K'^. 

A nonlinear complementarity problem (NCP) is composed 
of a number of nonlinear complementarity constraints || 6 l that 
take the form: 

a; > 0 ( 1 ) 

f(x) > 0 (2) 

= 0 (3) 

where a; G K"? and f —>■ K'?. 

A mixed linear complementarity problem (MLCP) is de¬ 
fined by the following constraints: 

Aa: -f Cy -I- g = 0 (4) 

Da; + By + h > 0 (5) 

y > 0 ( 6 ) 

y~^ {Cx + Dy + h) = 0 (7) 

Note that the x variables are unconstrained, while the y 
variables must be non-negative. If A is non-singular, the 


JOURNAL OF LTbX CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012 


2 


unconstrained variables can be computed as: 

x = -A^^{Cy + g) ( 8 ) 

Substituting x into equations |4j|^ yields the LCP (e,F): 

F = B - DA (9) 

e = h — DA^g (10) 

A solution (y, u) to this LCP obeys the relationship Fy+ e = 
i/; once one has y, x may be determined via Equation and 
the MCLP is solved. 


III. Background 

A. Coulomb friction 

Coulomb’s friction model provides relationships between 
the force applied along the contact normal and the frictional 
forces. Coulomb friction considers two cases, rolling/sticking 
and sliding. The former occurs when the velocity is zero in the 
tangent plane of the contact frame; conversely, sliding occurs 
when that velocity is non-zero. 

The magnitude of the friction force for a sliding contact 
modeled with Coulomb friction is given by the equation; 

ff=gcfn ( 11 ) 

where /„ is the magnitude of the force applied along the 
contact normal. The frictional force is applied directly opposite 
the direction of sliding (i.e., against the relative velocity in the 
tangent plane of the contact frame). 

The magnitude of the friction force for a rolling or stick¬ 
ing contact modeled with Coulomb friction is given by the 
equation: 

// < Mc/n (12) 

In the case of rolling/sticking friction, the friction force acts to 
resist motion in the tangent (e.g., in the case of a box resting 
on a slope); thus, ff may be strictly less than gcfn- If external 
forces become sufficiently large to overcome rolling/sticking 
friction forces, the rolling/sticking contact will transition to 
sliding. 

Many applications in robotics use the Coulomb friction 
model because it is relatively straightforward to compute— 
one can determine the frictional forces without integrating 
ordinary differential equations—and it is reasonably predic¬ 
tive. Nevertheless, Coulomb friction is somewhat expensive 
(computationally) to model: the rigid contact models of Stew¬ 
art and Trinkle ll23l and Anitescu and Potra HI can be 
solved in expected polynomial time in the n contact^ though 
exponential complexity may be exhibited in the worst case. 

B. Acceleration-level rigid body contact model with Coulomb 
and viscous friction 

We now describe the rigid contact model with Coulomb 
and viscous friction that uses only non-impulsive forces for 
consistency with the principle of constraints ifTOl . The multi 

* These models yield an order n copositive-plus LCP solvable by Lemke’s 
Algorithm 03 Each iteration of Lemke’s Algorithm requires an O(n^) 
matrix factorization update, and n iterations of the algorithm ai‘e expected (^. 


rigid body dynamics equation with contact and joint constraint 
forces is given below: 

M{t)v ^fif) + + mVfn + . . . (13) 

mtVfs + - ClifVfq - ... (14) 

S(f) Vs(f)-u - T(f) VT(t)r; (15) 


where M S 


is the system inertia matrix; v G 


pm 


Ujxm jjjg jnatrix of j bilateral 
S e K"''™, and T G 


is the system velocity; J G 
constraint equations; N G K 
are matrices of n wrenches applied along the contact normal, 
first contact tangent, and second contact tangent, respectively; 
Q G is a matrix of r generalized wrenches applied 

against the direction of sliding for r < n sliding contacts; fj G 
is the vector bilateral constraint force magnitudes; /„ G 
K" is a vector of contact normal force magnitudes; fg G 
and ft G are vectors of contact tangent force magnitudes 
applied at the k = n — r rolling/sticking contacts; fg G 
is a vector of contact tangent force magnitudes applied at the 
r sliding contacts; / is a vector of forces on the rigid body 
system (gravity, Coriolis and centrifugal forces, etc.); and /j,„ 
is a diagonal matrix of viscous friction coefficients. 

Equation specifies the Coulomb friction forces and 
Equation 15 specifies the viscous friction forces. The viscous 
model, where friction forces oppose the direction of motion, is 
commonly used in robotics (see, e.g., El) and is a simplifica¬ 
tion of the viscous drag term in fluid dynamics. Out of the n 
points of contact in the system, some may be rolling/sticking 
and the remainder will be sliding. Eor Coulomb friction, the 


first two terms of Equation 14 {Sk{tVfs + Tk(tV ft) 


relevant to the rolling/sticking contacts only (k specifies the 
indices of S and T that correspond to rolling/sticking contacts) 
and the last term (—Q(t)^/q) is relevant to only sliding 
contacts. 

1) Bilateral constraint equation: Bilateral constraints can 
be specified in the form (j){q) = 0, where q are the generalized 
coordinates of the system (joint constraints that are an explicit 
function of time are not considered here, though their inclusion 
would not change the results in this article). Such constraints 
can be differentiated once with respect to time to yield; 


Ji) = 0 (16) 

where J = ^. If we differentiate the constraints with respect 
to time once more, the bilateral joint constraints can be 
enforced using the equation: 

Jt) + jv = 0 (17) 

where J = ^3v. 

dq 

2) Contact normal constraints: We assume that there are 
n points of contact. The i* contact must satisfy the following 
linear complementarity condition that relates normal force and 
non-interpenetration: 

0 < fm -L rii^v -f rii^v > 0 (18) 

where nt and hi are column vectors taken from the rows 
of N and N, respectively. Here we adopt the notation a _L 6 
to denote the relationship a ■ b = 0. fm >0 requires that the 
contact force can only push bodies apart, nf^v -f nf^v > 0 
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requires that the bodies cannot be accelerating toward one 
another at the contact point after the contact forces are 
applied, and the /„. • [rij^v + v) = 0 constraint ensures 
that frictionless contact does no work. 

3) Sliding friction: If the velocity in the tangent plane at 
the contact point is non-zero, then the contact is sliding, 
and the Coulomb friction model specihes the magnitude of 
force to be applied. 

\U\=tt\fn.\ (19) 


4) Rolling/sticking friction: If the velocity at time t in 
the tangent plane at the z* contact point is zero, then the 
contact is rolling/sticking at time t and may either continue 
rolling/sticking or begin sliding, depending on the normal 
force and Coulomb friction coefficient. The nonlinear com¬ 
plementarity conditions are then expressed by the following 
equations (adapted from ll^ l: 


0 < uVn. - fl - 

fl -L > 0 

(20) 

+ 

II 

o 

/ 

(21) 

0 = BfriiVti + fti\i 

! 

(22) 


Let us now examine the constraints above. Equation 
constrains the frictional force to lie within the friction cone; if 
the tangential acceleration is non-zero, then the frictional force 
must lie on the edge of the friction cone. Equations and 22 


ensure that the frictional force opposes the tangent accelera¬ 
tion. 

5) Solvability of the model: Others {e.g., Q) have already 
shown that this rigid model may not possess a solution if 
there are any sliding contacts; such contact scenarios are 
known as inconsistent configurations. Nevertheless, we present 
this model because the contact model with Coulomb friction, 
which we present next and use to motivate the move to a 
velocity-level contact model, will build off of it. 


C. Solvable rigid body contact model with Coulomb and 
viscous friction 

The contact model of Stewart and Trinkle ll23l and Anitescu 
and Potra JT] provides a guaranteed solution to the problem of 
inconsistent configurations in contact models with Coulomb 
friction. This model is presented to show a velocity-level 
formulation, which allows the model to overcome the issue 
of inconsistent configurations. The no-slip model introduced 
in Section VI will also employ a velocity-level formulation 
to simulate contact without sliding in arbitrary configurations; 
Lynch and Mason showed that sliding with inhnite friction 
is possible for the acceleration-level model described in the 
previous section ns. 

We now describe this contact model—we consider only the 
aspect of the model that treats all contacts as inelastic impacts 
and do not consider extensions to collisional impacts with 
restitution. Eor simplicity of presentation, we do not linearize 
the friction cone, which yields a NCP rather than the LCP 
in lUmi]. 

The contact model uses a first-order approximation to the 
rigid body dynamics to resolve issues like Painleve’s Paradox 


(and other inconsistent contact configurations 121), for which 
no non-impulsive force solutions exist. The rigid body dynam¬ 
ics are given by: 

M.{t)lSv^lStf{t)+3{tf fj+... (23) 

ft 


where M, v, J, N, S, T, fn, fs, ft, and f are as defined 
in Section III-B and At is the change in time that realizes the 
first-order approximation. We now define v* = v + Av. 

1) Bilateral constraint equation: Because the bilateral joint 
constraints are now defined at the velocity level, the constraints 
are enforced using the equation: 


Jr;* = 0 (24) 

2) Contact normal constraints: The velocity-level con¬ 
straints on contact normal force and non-interpenetration are 
now defined as: 


0<fn^AnJv*>0 (25) 

3) Coulomb friction constraints: Coulomb friction is ef¬ 
fected more simply in this model than in the acceleration- 
level model: contacts can be treated identically whether they 
are initially sliding or sticking. The nonlinear complementarity 
conditions for the z* contact are: 

0 < u^fl - fl - fl T > 0 (26) 

0 = /r/„. < + /,. (27) 

0 = /r/„. vl + fu \l<+< (28) 

These equations are analogous to Equations p0] - |22] 

If the nonlinear complementarity conditions are converted 
to linear complementarity constraints by use of a linearized 
friction cone {i.e., a friction polygon), then a provably solvable 
copositive-plus LCP m results. However, Lemke’s Algo¬ 
rithm ifTSll is currently the only algorithm provably capable 
of solving copositive-plus LCPs. Lemke’s Algorithm can ex¬ 
hibit exponential complexity ua, though polynomial time is 
expected. 


IV. Contact model with purely viscous friction 


We now describe the contact model with purely viscous 
friction. We start from the multi rigid body dynamics equation 
at the acceleration level (Equations and [TS] ), which are 
reproduced below: 

M{t)v =f{t) + 3{tffj + -f ... 

S{t)^ ^yS{t)V — T(t)^ IIyT{t)V 


To this we add bilateral constraints (Equation |III-B 1 1 and 
the normal contact compressive force and non-interpenetration 


and complementarity constraints (Equation III-B2 1 , again re¬ 
produced below: 

Jz; -P jt; = 0 

0 < rii^v -P rii^v _L fn. >0 for z = 1,..., zz 
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Combining these equations yields the following MLCP; 


denoted by indices I and D, respectively, as follows: 


'M 

-JT 



i) 


'f 


o' 

J 

0 

0 


fj 

+ 

Jv 

= 

0 

N 

0 

0 


fn 


Nv_ 


.T 


fu>0 
7 > 0 
fnl = 0 


(29) 

(30) 

(31) 

(32) 


N = 


N/ 

Nd 


(43) 


Then the LCP vectors q = Nr;, 2 G K”, and u; G ffi" and 
LCP matrix Q = can be partitioned as follows: 


Qii 

Qdi 


Qid 

Qdd 



Zl 

+ 

qi 

_ 

Wj 


.^D_ 


. 9 d . 


Wd_ 


(44) 


where f* = — f + S~^fj,^Sv + T~^fiyTv and 7 = Nv + Nr;. 
As long as J has full row rank (we will describe how to ensure 
this condition in the next section), the mixed LCP can be 
converted to a conventional LCP (as described in Section 0 
using the following definitions: 


A = 


M 

J 


-JT 

0 



D = -C^ 
B = 0 


V 



y = fn 



h = 0 


(33) 

(34) 

(35) 

(36) 

(37) 

(38) 

(39) 

(40) 


Equations 1^ and 10 then yield the following standard LCP: 


F = NA^^N"^ (41) 

e = NA-^f* (42) 


The system inertia matrix is block diagonal (each block is 
invertible), so A is invertible if J has full row rank (if it is 
not—indicating that one or more constraints is redundant— 
a subset of J which has full row rank can be used to 
ensure that A is invertible). From 0 , a matrix of F’s 
form must be non-negative definite, i.e., either positive semi- 
definite (PSD) or positive definite (PD). Additionally, Baraff 
provided an algorithm that provably solved LCPs of the form 
(Gr,GHGT), where H G is a symmetric matrix, 

r G K"*, and G G Q. Finally, we note that LCPs 

with PSD/PD matrices are equivalent to convex quadratic 
programs which means that solving the LCP exhibits 
worst-case polynomial computational complexity. 


Given some matrix a G it is the case that = 

aN/, and therefore that Q^/ = aN/A^^N/^, Q/d = 
N/A^^N/^a^ (by symmetry), Qdd = aN/A^^N/^a^, 
and qo = aN/v. 

Lemma 1. Since rank('X.'Y) < min {rank{iX.), rank(Y)), the 
number of positive components of Zj can not be greater than 
rank{A). 

Proof Since the columns of XY have X multiplied by 
each column of Y, i.e., XY = [Xyi Xy 2 Xy„]. 

Columns in Y that are linearly dependent will thus produce 
columns in XY that are linearly dependent (with precisely the 
same coefficients). Thus, rank(XY) < rank(Y). Applying the 
same argument to the transposes produces 
rank(XY) < rank(X), thereby proving the claim. □ 

We now show that—in the case that the number of positive 
components of Zi is equal to the rank of A—no more positive 
force magnitudes are necessary to solve the LCP. 

Theorem 1. If {zj = a, wj = 0) is a solution to the 
LCP (g/,Q//), then {[zi~^ = ^d^ = 0^]^ ,w = 0) is 

a solution to the LCP {q, Q). 

Proof. For ([z/^ = aJ zd^ = 0^] ^ m = 0) to be a solu¬ 
tion to the LCP (q, Q), six conditions must be satisfied: 

1) z/ > 0 

2 ) wj > 0 

3 ) zj^wj = 0 

4) zd>0 

5) wd > 0 

6 ) zd^wd = 0 

Of these, ( 1 ), ( 4 ) and ( 6 ) are met trivially by the assumptions 
of the theorem. Since zd = 0, Qj/Z/ + QidZd + qi — 0, 
and thus wj — 0, thus satisfying ( 2 ) and ( 3 ) Also due to 
Zd = 0 , it suffices to show for (5) that QdiZi + qo > 0 . 
From above, the left hand side of this equation is equivalent to 
a (N/A-^N/'^a -t- N/d), or awj, which itself is equivalent 
to aO. Thus, Wo = 0. □ 


V. Reducing expected time complexity 

FROM O(n^) TO Ofnf + mfn) 

A system with m degrees-of-freedom requires no more than 
m positive force magnitudes applied along the contact normals 
to satisfy the constraints for the contact models with purely 
viscous friction and without slip (the latter model will be 
presented in Section [VU l. We now prove this statement. 

Assume we permute and partition the rows of N into 
r linearly independent and n — r linearly dependent rows. 


4} Algorithm: We use the theorem above to make a mi¬ 
nor modification to the Principal Pivot Method 1 121 m 
(PPM), which solves LCPs with P-matrices (complex square 
matrices with fully non-negative principal minors ifTSll that 
includes positive semi-definite matrices as a proper subset). 
The resulting algorithm limits the size of matrix solves and 
multiplications. 

The PPM uses a set /? with maximum cardinality n for a 
LCP of order n. Of a pair of LCP variables, {zi,Wi), exactly 
one will be in /3; we say that the other belongs to f3. If a 
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variable belongs to 13, we say that the variable is a basic 
variable', otherwise, it is a non-basic variable. Using this set, 
partition the LCP matrices and vectors as shown below: 


Wp 


App App 


'zp 

_l_ 




A_ A_ 

L^/3/3 ^/3/3j 




.9/3. 


Segregating the basic and non-basic variables on different sides 
yields: 


Wp 


^PP ^pp^pp ^pp-^pp 


^p 



_A— _ A_ 




qp - 

If we set the values of the basic variables to zero, then solving 
for the values of the non-basic variables and uig entails 
only block inversion of A. 

The unmodihed PPM I operates in the following manner: 

(1) Find an index i of a basic variable Xi (where Xi is either Wi 
or Zi, depending which of the two is basic) such that Xi < 0 ; 

( 2 ) swap the variables between basic and non-basic sets for 
index i (e.g., if Wi is basic and is non-basic, make Wi non- 
basic and Zi basic); ( 3 ) determine new values of 2 ; and w; ( 4 ) 
repeat (i)-( 3 ) until no basic variable has a negative value. 

PPM I requires few modihcations, which are provided in 
Algorithm Iri First, the full matrix N • is never 

constructed (such construction would require 0{n^) time). 
Instead, Line [T] of the algorithm constructs a maximum mxm 
system; thus, that operation requires only 0{m?) operations. 
Similarly, Lines 13-14 also leverage Theorem [T] in order to 
compute w''' and efficiently (though these operations do 
not affect the asymptotic time complexity). Assuming that 
the number of iterations for a pivoting algorithm is 0{n) 
in the size of the input]^ and that each iteration requires at 
most two pivot operations (each rank -1 update operation to 
a matrix factorization will exhibit time complexity 0{m?)), 
the asymptotic complexity of the modified PPM I algorithm 
is 0{m^ +'m?n). The termination conditions for the algorithm 
are not affected by our modihcations. 

VI. No-SLIP CONTACT MODEL 

A contact model without slip requires a velocity-level 
contact model in accordance with Lynch and Mason’s hnding 
that sliding can occur with inhnite friction at the acceleration 
level El. The no-slip friction contact model uses the hrst- 
order approximation and builds on Equations and 

by dictating that the tangential velocity at each contact must 
be zero at v{t-\- At): 


Forming these hve equations into a MLCP and using the 
variable dehnitions in Section [H] yields: 


A = 


M 

J 

S 

T 


-JT 

0 

0 

0 


■-NT‘ 

c= ° 

0 

D = -C^ 

B = 0 


_ST _-pT 

0 0 

0 0 

0 0 


(47) 


(48) 

(49) 

(50) 

(51) 


v(t + At) 

0 

0 


y=fn 


g = 


-Mv{t) 

0 

0 

0 


h = 0 


(52) 

(53) 

(54) 

(55) 


The matrix A may be singular, which would prevent us from 
converting the MLCP to a standard LCP. However, if J, S, 
and T have full row rank or we identify the largest row 
blocks of those matrices such that full row rank is attained, 
A is invertible without affecting the solution to the MLCP. 
Algorithm 1^ performs exactly this task. 

The singularity check on Lines]? 14 and 19 of Algorithm]^ 
is best performed using a Cholesky factorization; if the fac¬ 
torization is successful, the matrix is non-singular. Given that 
M is non-singular (it is symmetric and positive definite), the 
maximum size of X in Algorithm is m x m; if X were 
larger, it would be singular (see Lemma [T]). 

Given this information, the time complexity of Algorithm]^ 
is dominated by Lines IZlEl and As X changes by at 
most one row and one column per Cholesky factorization, 
singularity can be checked by 0{mf) updates to an initial 
0{rnf) Cholesky factorization. The overall time complexity 
is 0{m^ -f nmf). 


A. Resulting systems 

Using Equations and [T0] l, the LCP matrix F and vector 
e are equivalent to: 


Sv{t + At)=0 (45) 

Tv{t + At)=0 (46) 


F = NX-^N't (56) 

e = NX~^Mv{t) (57) 


^Regardless of the pathological problem devised by Klee and Mintv illl . 
experience with the Simplex Algorithm on thousands of practical problems 
shows that it requires fewer than Sn iterations and the expected time 
complexity for the Simplex Algorithm is polynomial I18II20I . We are unaware 
of research that shows these results are also applicable to pivoting methods 
for LCPs, though Cottle et al. claim 0(n) expected iterations O. 


As in Section |IV| F must be symmetric and positive-semi- 
dehnite and—as noted in Section jiy]— Baraff’s algorithm ||2| 
guarantees that a solution to this LCP exists. 

The Sv{t + At) = Tv{t + At) = 0 constraints (Equa¬ 
tions 1^ and |46]l and solvability of the LCP contrast with the 
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Algorithm 1 {z,w} = LCP{N,M, f*) Solves a frictionless contact model using a modification of the Principal Pivoting 
Method I Algorithm. 

1 

n 3— rows(N) 


2 

q^N-r 


3 

i 4— arg miUj qi 

[> Check for trivial solution 

4 

if gi > 0 then 


5 

return { 0 , < 7 } 


6 

end if 


7 

B 3 - {z} 

> Establish initial nonbasic indices 

8 

B ^ {1,... ,i - l,i + 1,... ,n} 

> Establish initial basic indices 

9 

while true do 


10 

A ^ Nb • M-i • 


11 



12 

^ A-i • -b 

\> Solve for non-basic components of 2 : 

13 

at ^ M-i • Ng'^^t + f* 


14 

iijt ^ N • at 


15 

i 4— arg minj wj 

> Find the index for moving into the non-basic set (if any) 

16 

if wj > 0 then 


17 

j <— arg miuj z^ 

> No index to move into the non-basic set; look whether there is an index to 

18 


move into the basic set 

19 

if z] < 0 then 


20 

k ^ B(j) 


21 

B^BU{k} 

[> Move index k into the basic set 

22 

B^B-{k} 


23 

continue 


24 

else 


25 

2 ; 3 — 0 


26 

2g^2;t 


27 

uj 3 — 0 


28 

lUg 3— mt 


29 

return {z, to} 


30 

end if 


31 

else 


32 

B^BU{i} 

> Move index i into the non-basic set 

33 

B^B-{i} 


34 

j 4 — arg miuj 2 ;t 

> Look whether there is an index to move into the basic set 

35 

if 2 ] < 0 then 


36 

k ^ B{j) 


37 

B^BU{k} 

[> Move index k into the basic set 

38 

B^B-{k} 


39 

end if 


40 

end if 


41 

end while 



finding of Lynch and Mason m, who showed that sliding 
with infinite friction is possible. The admittance of impulsive 
forces has resolved this “paradox” analogously to the manner 
in which contact models like ll2^ fTll resolved Painleve’s Para¬ 
dox HU and other inconsistent contact configurations Il22l . 

VII. Experiments 

We tested the contact models using two common con¬ 
tact scenarios in robotics, grasping and locomotion, in or¬ 
der to assess speed and numerical stability. These experi¬ 
ments can be reproduced using the experimental setup de¬ 
scribed at https://github.com/PositronicsLab/ 
no-siip-and-viscous-experiments 


A. Grasping experiment 

We used RPIsim (https://code.google.eom/p/ 
rpi-matlab-simulator) to simulate a force-closure 
grasping scenario (depicted in Figure [^1 on a Macbook Air 
with 1.8 GHz Intel Core i5 CPU. Twelve contact points were 
generated between each pair of boxe^ yielding 36 contact 
points total. For the contact model with Coulomb friction 
(the Stewart-Trinkle model El), a friction “pyramid” (four 

^The equal size boxes contacting in the manner in Figure ^ yields 
degenerate contact normals at the box comers; the RPI simulator treats this 
problem by duplicating each contact with all three possible directions for the 
contact normal. 
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Algorithm 2 Find-Indices(M, J, S, T), determines the row 
indices { J, S, and T) of J, S, and T such that the equality 
matrix A (Equation is non-singular. 

1 : 

2 : 

3: r^0 

4: for J = 1, ... r do > r is the number of bilateral 

constraint equations 
5: 

6: Set X ^ J}, 

7: if not singular then 

8: J ^ J* 

9: end if 

10 : end for 

11 : for i = 1,... ,n do > n is the number of contacts 

12 : 5*4-5U{2} 

13: Set X ^ \3'lj Sj. T^] 

14: if X^M^^X not singular then 

15: 5^5* 

16 : end if 

17: T*^TU{t} 

18: Set X ^ [jT sj t:^.] 

19: if X^M~^X not singular then 

20: T ^T* 

21 : end if 

22 : end for 

23: return {J,S,T} 


sided approximation to the friction cone) was used, yielding 
six LCP variables per contact (i.e., 216 variables total). The 
RPI simulator allowed us to substitute the Stewart-Trinkle 
model with the no-slip friction model readily, which resulted 
in only 36 LCP variables. The simulation was run using 
a step size of 0.01 for ten iterations of one second 
of simulated time); the grasped objects would tend to fall 
from the gripper after ten iterations only when using Stewart- 
Trinkle (due to numerical issues with Lemke’s Algorithm, to be 
discussed below). Lemke’s Algorithm was implemented using 
LEMKE 19). 

The Pq matrix resulting from the no-slip friction model 
allowed us to employ the modified PPM solver and MATLAB’s 
quadprog solver (with the active-set algorithm) to solve the 
LCP. We used Lemke’s Algorithm ifT^ . employing Tikhonov 
regularization ID as necessary, to solve the Stewart-Trinkle 
model. No low-rank updates were used in our implementation 
of Algorithmic 


Contact model 

Running time (mean di a) 

Stewart-Trinkle (Lemke’s Algorithm) 

flc = 100.0, /2i; = 0.0 

No-slip (active-set QP solver) 

No-slip (modified PPM) 

10.9681s ± 2.1812s 

1.9892s ± 0.2640s 
1.6680s ± 0.3669s 


TABLE I 

Mean running times eor the grasping experiment. Ten trials 

WERE RUN FOR EACH METHOD. TIMINGS INCLUDE ALL ASPECTS OE THE 
SIMULATION (INCLUDING COLLISION DETECTION). 


This experiment yielded several findings. As expected. 



Fig. 1. A depiction of the grasping experiment described in Section [Vll-A| 
The two red boxes act as grippers and push inward. Gravity pushes downward. 
Given sufficient friction (jt = oo), the grippers should ideally keep the blue 
boxes grasped using force closure. 


reducing the LCP variables by a factor of five (216 variables 
to 36 variables) results in much faster solutions (451-558% 
faster mean, depending on the solver). Lewer variables also 
results in less rounding error; the no-slip approach was able to 
model the grasping scenario reliably for at least 100 iterations 
(again, compared to around ten iterations for Stewart-Trinkle). 
Of the twelve contacts per pair of boxes, it was only necessary 
to apply forces to two contacts, which the modified PPM 
method was able to exploit; it ran nearly 20 % faster than 
the quadprog algorithm (mean and maximum numbers of 
pivot operations were observed to be 5.5 and 7, respectively). 
This performance differential is considerable given that our 
modified PPM algorithm was not implemented as a MEX file 
and that our implementation does not use low-rank updates to 
maximize performance). 

B. Locomotion experiments 

We used the Moby simulator (https : //github . com/ 
PositronicsLab/Moby) to simulate a quadrupedal robot 
walking on a terrain map (see Ligure ^ over ten seconds. 
An event-driven method (see i) for a description of this 
paradigm) is used to simulate the system instead of the time¬ 
stepping approaches used in ODE, Bullet, and RPIsim; 
popular implementations of this approach are susceptible 
to energy gain when correcting interpenetration ll2Tll . which 
destabilizes our robot in the process. 

The integration method used for the no-slip model experi¬ 
ment is symplectic Euler (Stormer-Verlet) with a step size of 
0.001, while fourth-order Runge-Kutta integration was used 
for the viscous model experiment with identical step size. 
Our approach using the former integrator allows only the 
no-slip model to be activated, while our approach using the 
latter integrator permits the acceleration-level viscous model 
to be used for sustained contacts. When impacts occur (e.g., 
on initial foot/ground contact) in the latter approach, the 
simulation uses an inelastic impact model 13 with purely 
viscous friction. 

Locomotion experiments were run on an Intel Xeon 
2.27GHz desktop computer. All aspects of the simulation, 
including forward dynamics computations, collision detection. 
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and controls (which performs dynamics calculations) were 
accounted for in results; single LCP solves generally run 
too quickly to obtain timings for just that operation. We 
point out that the simulations run considerably slower than 
similar systems modeled using, e.g., ODE, but the goals of 
the two simulators. ODE uses approximate solves and permits 
interpenetration. Moby aims to provide a verifiable simulator, 
i.e., one that adheres closely to the rigid body dynamics models 
(which means that interpenetration is impermissible). 

Results for the no-slip experiment are provided in Ta- 
which shows a speedup of nearly 28%. 


ble Vll-B 


The 

minimum and maximum number of contact points generated 
in the experiment is 1 and 30, respectively; the mean number 
of contact points is 6, and the standard deviation is 8. Thus, 
simulations with greater numbers of contacts could expect 
greater performance differentials. 

Results for the viscous experiment are provided in Ta¬ 
ble VITB which shows a speedup of over 37%. We note 
that the viscous friction experiments required significantly 
longer to run than the no-slip experiments. We hypothesize 
that this disparity is due to the behavior of the simulation 
when applying the viscous model, which tends to produce 
rapid movements upon contact. Those rapid movements, which 
appear due to some sensitivity in the underlying ordinary dif¬ 
ferential equations, slow the simulator’s continuous collision 
detection system (see m for a description of that system). 



Fig. 2. A depiction of the quadrupedal robot walking on a tenain map in 
the locomotion experiment, as described in Section [Vll-B| The no-slip and 
viscous friction models were both assessed. 


Contact model 

Running time 

Drumwright-Shell 

fic = 100.0, fLv = 0.0 

No-slip (modified PPM) 

416.284s 

301.796s 


lABLH II 


Times required to simulate the quadruped locomotion 
SCENARIO DESCRIBED IN SECTIOn IVII-BI UNDER A NO-SLIP CONTACT 
MODEL. Timings include all aspects of the simulation 

(INCLUDING collision DETECTION). 


VIII. Conclusion 

We presented an algorithm for rapidly computing two rigid 
contact models without Coulomb friction that have proven 
useful in certain modeling and simulation applications for 


Contact model 

Running time 

Viscous (Lemke’s Algorithm) 
/ic = 0.0, fiv = 0.1 

Viscous (modified PPM) 

2936.36s 

2139.73s 


TABLE 111 

Times required to simulate the quadruped locomotion 
scenario described in Section IviI-BI under a purely viscous 
friction model. Timings include all aspects of the simulation 
(including collision detection). 


robotics. We showed how these models exhibit both asymp¬ 
totic computational complexity and signihcant running time 
advantages over rigid models with Coulomb friction. While we 
do not expect these special-case models to replace rigid models 
with Coulomb friction, the former serve as computationally 
efficient alternatives as applications allow. 
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